#include "stdafx.h"
//#include "RothC_Input.h"
//
// Par_Maps RothC::par_rothc_crop_pool_distri;
// Par_Maps RothC::par_rothc_fert_pool_distri;
// Par_Maps RothC::par_fert_props;
// Par_Maps RothC::par_rothc_EWR_franko;
// Par_Seed_Harvest RothC::par_default_date_seeding_harvest;
//RothC::RothC(RothC_Input &driver)
//{
//    double best_fit;
//
//    //Pools	Indices
//    dpm=0;rpm1=1;rpm2=2;hum=3;bio=4;iom=5;
//	
//    // Response function indices
//    respT=0;respWC=1;respVeg=2;
//    // non_CO2_fraction
//    if((driver.plotpar["clay"]<0)||(driver.plotpar["clay"]>100)) err::err2("RothC::RothC: clay out of Range");
//    	non_CO2_fraction=1/(1+1.67*(1.85+1.6*exp(-0.0786*driver.plotpar["clay"])));
//
//    // set default pool amounts
//    init_est_pool_stocks ( 0,driver);
//
//    // call equilibrium run
//   best_fit= call_equlibrium_run(driver);
//   time=driver.finaltime;
//  
//}
//
//int RothC::init_est_pool_stocks( int iom_mode, RothC_Input &driver)
//{
//    double initSOC;
//    if ((driver.plotpar["SOC"]<0)||(driver.plotpar["SOC"]>100)) err::err2("RothC::init_est_pool_stocks: initSOC out of Range");
//    if((driver.plotpar["bd"]<0)||(driver.plotpar["bd"]>3))err::err2("RothC::init_est_pool_stocks: bd out of range.");
//    //conversion from [%] SOC into [t/ha]
//    initSOC=driver.plotpar["SOC"]/100*driver.plotpar["bd"]*100*100*0.3;
//    switch (iom_mode)
//	{
//		//Falloonequation
//		case 0:
//                   
//                    Cpool [0][iom]=0.049*pow((double)initSOC,(double)1.139);
//			break;
//		//Körschens equation
//                        //iom=0.04*bd[g/cm3]*clay[%]*soil_sample_depth;
//		case 1:
//                    if((initSOC<0)||(initSOC>999))err::err2("RothC::init_est_pool_stocks: initSOC out of range.");
//                   
//                    if((driver.plotpar["clay"]<0)||(driver.plotpar["clay"]>100))err::err2("RothC::init_est_pool_stocks: clay out of range.");
//                    if((driver.plotpar["sample_depth"]<=0)||(driver.plotpar["sample_depth"]>50))err::err2("RothC::init_est_pool_stocks: sample depth out of range.");
//
//                    Cpool[0][iom]=0.04*driver.plotpar["clay"]*driver.plotpar["bd"]*driver.plotpar["sample_depth"];
//
//                    if((Cpool[0][iom]<0)||(Cpool[0][iom]>=initSOC)) err::err2("RothC::init_est_pool_stocks: IOM out of Range.");
//			break;
//                default: err::err2("RothC::init_est_pool_stocks: wrong identifier that determines the way of iom estimation.");
//                    break;
//	}
//            /*rate..DPM:*/ Cpool[0][dpm]=0.01*(initSOC-Cpool[0][iom]);
//            /*rate..RPM1:*/ Cpool[0][rpm1]=0.06*(initSOC-Cpool[0][iom]);
//            /*rate..RPM1:*/ Cpool[0][rpm2]=0.06*(initSOC-Cpool[0][iom]);
//            /*rate..BIO:*/ Cpool[0][bio]=0.02*(initSOC-Cpool[0][iom]);
//            /*rate..HUM:*/ Cpool[0][hum]=0.85*(initSOC-Cpool[0][iom]);
//      return 0;
//}
//
//double RothC::call_equlibrium_run( RothC_Input &driver)
//{
//	unsigned int run;
//	double ax=0,bx=1,cx=15,tol=0.000001,corg=0;
//	double *xmin;
//	double input;
//	map<const string, vector <double> >::iterator it;
//	vector <double> croptype;
//
//	xmin=&corg;
//	timestep=365.25/12.0;
//
//	if (driver.tempdisc->size()!=12)err::err2("Driving date of equilibrium run has less or more than 12 data entries.");
//
//        //Tresponse
//        if(driver.met["T"].size()!=12)err::err2("RothC::call_equilibrium_run:Driving date of equilibrium run has less or more than 12 data entries");
//        it=driver.inp.find("Tresponse");
//	if (it!=driver.inp.end())  it->second.clear();
//        for(run=0;run<driver.inp["T"].size();run++) driver.inp["Tresponse"].push_back(get_Tresponse(driver.inp["T"][run]));
//
//        //WCResponse
//        get_wc(driver);
//        if(driver.inp["WC"].size()!=12)err::err2("RothC::call_equilibrium_run:Driving date of equilibrium run has less or more than 12 data entries");
//        it=driver.inp.find("WCresponse");
//	if (it!=driver.inp.end())  it->second.clear();
//        for(run=0;run<driver.inp["WCrresponse"].size();run++) driver.inp["WCresponse"].push_back(get_WCresponse(driver.inp["WC"][run],get_maxwc(driver.plotpar["clay"])));
//
//        //Soilcover response
//        it=driver.inp.find("crop_for_equi_run");
//	if (it!=driver.inp.end())  it->second.clear();
//	par_default_date_seeding_harvest.Annual_time_line((int) driver.plotpar["cropid_for_equi_run"],NOCROP,croptype);
//        driver.inp["crop_for_equi_run"]=croptype;
//        it=driver.inp.find("SOCOresponse");
//	if (it!=driver.inp.end())  it->second.clear();
//        for(run=0;run<driver.inp["crop_for_equi_run"].size();run++) driver.inp["SOCOresponse"].push_back(get_SOCOresponse(driver.inp["crop_for_equi_run"][run]));
//
//	input=brent(ax,  bx,  cx, tol, xmin,driver);
//	return input;
//}
//
//double RothC::call_RothC_during_equilibrium_run (RothC_Input &driver, double iCinput)
//{
//	int y=0,month=0;
//	int fertid,cropid;
//	double corg_sum_new,opti;
//
//	fertid=(unsigned int)driver.plotpar["fertid_for_equi_run"];
//	cropid=(unsigned int)driver.plotpar["cropid_for_equi_run"];
//      // partitioning of the Cinput between pools
//        Cinput[dpm]=iCinput*par_rothc_crop_pool_distri.get_val("dpm",cropid)+
//                par_rothc_fert_pool_distri.get_val("dpm",fertid);
//        Cinput[rpm1]=0.5*iCinput*par_rothc_crop_pool_distri.get_val("rpm",cropid)+
//                0.5*par_rothc_fert_pool_distri.get_val("rpm",fertid);
//        Cinput[rpm2]=0.5*iCinput*par_rothc_crop_pool_distri.get_val("rpm",cropid)+
//                0.5*par_rothc_fert_pool_distri.get_val("rpm",fertid);
//        Cinput[hum]=iCinput*par_rothc_crop_pool_distri.get_val("hum",cropid)+
//                par_rothc_fert_pool_distri.get_val("hum",fertid);
//        for(y=0;y<=3;y++) Cinput[y]=Cinput[y]/12.0;
//        do
//        {
//
//            response[respT]=driver.inp["Tresponse"][month];//temp response
//            response[respWC]=driver.inp["WCresponse"][month];//wfps response
//            response[respVeg]=driver.inp["SOCOresponse"][month];// soilcoverage response
//            y++;month=(month<11)*(1+month);
//            corg_sum_new=model();
//	}
//	while (y<=EQUICONST);
//	opti=(corg_sum_new-driver.plotpar["SOC"])*(corg_sum_new-driver.plotpar["SOC"]); //Berechnung der Modellanpassung
//	//printf ( "%i opti %f %f %f %f\n",x,opti,corg_sum_new,fi,c_input);
//	return opti;
//}
//
//double RothC::call_RothC(RothC_Input &driver)
//	{
//
//            double Cinput_crop;
//            double Cinput_fert;
//            double corg_sum_new;
//            unsigned int endrow,row0;
//            int y;
//            timestamp actual_time;
//            // check if itime is in inp
//            endrow=Basics_vector::find_position_of_object(driver.inptime,driver.finaltime); //actual time
//            if(endrow==0) err::err2("RothC::callRothC: itime is equal to adjacent time");
//            if(endrow==driver.inp["year"].size())err::err2("RothC::timestamp matches no step in driver time series");
//            y= row0=Basics_vector::find_position_of_object(driver.inptime,time); //previous time
//            if(row0>=endrow) err::err2("RothC::adjacent time does not occur in driver time series");
//           
//            //check for wc;
//            if(driver.inp.find("WC")==driver.inp.end()) get_wc(driver);
//
//            do
//            {
//                y=y+1;
//                actual_time=driver.inptime[y];
//                if (actual_time>driver.finaltime) err::err2("RothC::call_RothC: Final time is not defined in input time series.");
//                timestep=time.get_dist(actual_time);
//            //set Cinput zero
//                set_Cinput(0.0);
//                Cinput_crop=0;
//                Cinput_fert=0;
//            // get input from EWR
//                    switch(EWR_MODE)
//                    {
//                            case 0:Cinput_crop=get_EWR_franko((int)driver.inp["cropid"][y],driver.inp["cropyield"][y]);break;
////                            case 1:;break;
////                            case 2: /*klimanek_part_crop(x,y,month)*/;break;
////                            case 3: /*harvest_done=crop_growth2(x,y,month,harvest_done)*/;break;
//                            default:err::err2("RothC::call_RothC: EWR_MODE is ill defined");break;
//                    }
//              // partition EWR Input into RothC Pools
//                    Fill_Cinput_with_EWR((int)driver.inp["cropid"][y],Cinput_crop);
//             // get input from fertilizer
//                    Cinput_fert=get_C_from_fert((int)driver.inp["fertid"][y],driver.inp["fertamount"][y]);
//             //partition fertilizer into RothC pools
//                    Fill_Cinput_with_fert((int)driver.inp["fertid"][y],Cinput_fert);
//
//
//                     response[respT]=get_Tresponse(driver.inp["T"][y]);//temp response
//                     response[respWC]=get_WCresponse(driver.inp["WC"][y],get_maxwc(driver.plotpar["clay"]));//wfps response
//                     response[respVeg]=get_SOCOresponse((int)driver.inp["cropid"][y]);// soilcoverage response
//            
//            // start carbon model
//                    corg_sum_new=model();
//                
//
//                    time=actual_time;
//            }
//            while (time<driver.finaltime);
//           
//            return corg_sum_new;
//	}
//
//
//double RothC::brent(double ax, double bx, double cx, double tol,double *xmin,RothC_Input &driver)
////Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is
////between ax and cx, and f(bx) is less than both f(ax) and f(cx)), this routine isolates
////the minimum to a fractional precision of about tol using Brent's method. The abscissa of
////the minimum is returned as xmin, and the minimum function value is returned as brent, the
////returned function value.
//{
//int iter;
//int ITMAX=100;
//double CGOLD=0.3819660,ZEPS=0.01;
//double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
//double e=0.0;					//This will be the distance moved on
//									//the step before last.
//a=(ax < cx ? ax : cx);				//a and b must be in ascending order,
//b=(ax > cx ? ax : cx);		//input abscissas need not be.
//x=w=v=bx;				//Initializations...
////printf("hier bin ich noch\n");
//fw=fv=fx=call_RothC_during_equilibrium_run(driver,x);
////printf("hier nicht mehr");
//for (iter=1;iter<=ITMAX;iter++)//
//{// Main program loop.
//	//printf(" %i",iter);
//	xm=0.5*(a+b);
//	tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
//	//printf (" %i %f %f %f\n",iter, pool[0][2],pool[0][3],res_file[site][2]);
//	if (fabs(x-xm) <= (tol2-0.5*(b-a)))   //Test if the task is matched
//	{ //Test for done here.
//		*xmin=x;
//	//	res_file[site][2]=fx;
//	//	res_file[site][3]=constant[site][25];
//
//		return x;//eigentlich return fx
//	}
//	if (fabs(e) > tol1)
//	{ //Construct a trial parabolic t.
//		r=(x-w)*(fx-fv);
//		q=(x-v)*(fx-fw);
//		p=(x-v)*q-(x-w)*r;
//		q=2.0*(q-r);
//		if (q > 0.0) p = -p;
//		q=fabs(q);
//		etemp=e;
//		e=d;
//		if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
//		d=CGOLD*(e=(x >= xm ? a-x : b-x));
//		//The above conditions determine the acceptability of the parabolic t. Here we
//		//take the golden section step into the larger of the two segments.
//		else
//		{
//			d=p/q; //Take the parabolic step.
//			u=x+d;
//			if (u-a < tol2 || b-u < tol2)
//			d=((xm-x)>=0)*fabs(tol1)+((xm-x)<0)*-fabs(tol1);
//		}
//	}
//	else
//	{
//		d=CGOLD*(e=(x >= xm ? a-x : b-x));
//	}
//	u=(fabs(d) >= tol1 ? x+d : x+(((d)>=0)*fabs(tol1)+((d)<0)*-fabs(tol1)));
//	fu=call_RothC_during_equilibrium_run(driver,u);
//
//	//This is the one function evaluation per iteration.
//	if (fu <= fx)
//	{					//Now decide what to do with our funcif
//		if(u >= x) a=x; else b=x;			//tion evaluation.
//		v=w;w=x;x=u;	//SHFTHousekeeping follows:
//		fv=fw;fw=fx;fx=fu;//SHFT(fv,fw,fx,fu);
//	}
//	else
//	{
//		if (u < x) a=u; else b=u;
//		if (fu <= fw || w == x)
//		{
//			v=w;
//			w=u;
//			fv=fw;
//			fw=fu;
//		}
//		else if (fu <= fv || v == x || v == w)
//		{
//			v=u;
//			fv=fu;
//		}
//	} //Done with housekeeping. Back for
//	//printf("\n");
//} //another iteration.
//err::err2("Too many iterations in brent");
//return 0;
//}
//
//double RothC::model()
//   {
//       double degrad_sum=0, corg_sum_new=0;
//	double rrate[10],degrad[10];
//	int z;
//	// added Carbon
//	Cpool[0][dpm]=Cpool[0][dpm]+Cinput[dpm];//DPM
//	Cpool[0][rpm1]=Cpool[0][rpm1]+Cinput[rpm1];//RPM
//         Cpool[0][rpm2]=Cpool[0][rpm2]+Cinput[rpm2];//RPM
//	Cpool[0][hum]=Cpool[0][hum]+Cinput[hum];//HUM
//
//	/*rate..DPM:*/ rrate[dpm]=1-exp(-kC[dpm]*response[respT]*response[respWC]*response[respVeg]*timestep/365.25);
//	/*rate..RPM:*/ rrate[rpm1]=1-exp(-kC[rpm1]*response[respT]*response[respWC]*response[respVeg]*timestep/365.25);
//				   rrate[rpm2]=1-exp(-kC[rpm2]*response[respT]*response[respWC]*response[respVeg]*timestep/365.25);
//	/*rate..BIO:*/ rrate[bio]=1-exp(-kC[bio]*response[respT]*response[respWC]*response[respVeg]*timestep/365.25);
//	/*rate..HUM:*/ rrate[hum]=1-exp(-kC[hum]*response[respT]*response[respWC]*response[respVeg]*timestep/365.25);
//	/*rate..IOM:*/
//	// Degradation:
//		for (z=0;z<=4;z++)
//		{
//			Cpool[1][z]=Cpool[0][z]*(1-rrate[z]);degrad[z]=Cpool[0][z]*rrate[z];
//			degrad_sum=degrad_sum+degrad[z];
//		}
//
//	 // adding all the stuff to bio and hum
//		Cpool[1][bio]=Cpool[1][bio]+0.46*non_CO2_fraction*degrad_sum;
//		Cpool[1][hum]=Cpool[1][hum]+0.54*non_CO2_fraction*degrad_sum;
//	 // pool[1]-->pool[0]
//		for (z=0;z<=4;z++)
//		{
//			Cpool[0][z]=Cpool[1][z];corg_sum_new=Cpool[0][z]+corg_sum_new;
//		}
//		corg_sum_new=corg_sum_new+Cpool[0][4];
//	return corg_sum_new;
//   }
//
//double RothC::get_Tresponse(double temp)
// {
//     double alfa;
//// funktionsimplementierung
//		if (temp<-18)
//		{
//			alfa=0;
//		}
//		else
//		{
//			alfa= 47.9/(1+exp(106/(temp+18.3)));
//		}
//		//printf("%i month %i _sw %i %f alf %f\n",x,month,sw_,equi_inf[x][month][2],alfa);
//		return alfa;
// }
//
//double RothC::get_WCresponse(double wc, double maxwc)
// {
//    // funktionsimplementierung
//    double beta;
//    beta=(wc<0.444*maxwc)*1+((wc>=0.444*maxwc)&&(wc<maxwc))*(0.2+(1-0.2)*((maxwc-wc)/(maxwc-0.444*maxwc)))+(wc>=maxwc)*0.2;
//    //printf("%i month %i _sw %i %f beta %f\n",x,month,sw_,TSMD,beta);
//    return beta;
// }
//
//double RothC::get_maxwc (double clay)
//{
//    return (20+1.3*clay-0.01*pow(clay,2));
//}
//
//double RothC::get_SOCOresponse (double cropid)
//{
//    return((cropid==NOCROP)*1+(cropid!=NOCROP)*0.6);
//}
//int RothC::get_wc(RothC_Input &driver)
// {
// // WC Berechnen
//        map <const string,vector<double> >::iterator it;
//        it=driver.inp.find("WC");
//        if(it!=driver.inp.end()) it->second.clear();
//        if(driver.inp["P"].size()!=driver.inp["ETP"].size()) err::err2("RothC::get_wc:different number of P and T entries");
//        Basics_weather::get_bucket_storage(driver.inp,"P","ETP","WC",0,get_maxwc(driver.plotpar["clay"]));
//        return 0;
// }
//
//int RothC::get_bucket_storage(map<const string,vector<double> > &weather,string source,string sink, string target,double minvol, double maxvol)
//	{
//		double wc_init;
//		unsigned int row;
//		//Tests:
//		if((weather[source].size()<1)||(weather[sink].size()<1)||(weather[target].size()!=0))err::err2("Basics_weather::bucket_content: input map containing sources,sinks,target has strange size");
//		wc_init=maxvol;
//		for (row=0;row<weather[source].size();row++)
//		{
//			if((weather[source][row]<0)||(weather[sink][row]<0))err::err2("Basics_weather::bucket_content: Either source or sink value is negative");
//			if (row==0){
//				weather[target].push_back(Basics::inter(minvol,wc_init+(weather[source][row]-weather[sink][row]),maxvol));
//			}
//			else{
//				weather[target].push_back(Basics::inter(minvol,weather[target][row-1]+(weather[source][row]-weather[sink][row]),maxvol));
//			}
//		}
//		return 0;
//	}
//
//double RothC::get_EWR_franko(int cropid,double yield)
//  {
//      double a,b;
//      a=par_rothc_EWR_franko.get_val("a",cropid);
//      b=par_rothc_EWR_franko.get_val("b",cropid);
//      return a+b*yield; // [t/ha]
//
//  }
//double RothC:: get_C_from_fert(int fertid,double amount)
// {
//    return amount*par_fert_props.get_val("Dry_matter_content",fertid)*par_fert_props.get_val("C_content",fertid);
// }
//int RothC::set_Cinput(double a)
// {
//     int row;
//     for(row=0;row<NPOOLS;row++)
//     {
//         Cinput[row]=a;
//     }
//     return 0;
// }
//int RothC::Fill_Cinput_with_fert(int fertid,double Fert_C)
// {
//        Cinput[dpm]= Cinput[dpm]+par_rothc_fert_pool_distri.get_val("DPM",fertid)*Fert_C;
//        Cinput[rpm1]=Cinput[rpm1]+par_rothc_fert_pool_distri.get_val("RPM",fertid)*Fert_C*0.5;
//        Cinput[rpm2]=Cinput[rpm2]+par_rothc_fert_pool_distri.get_val("RPM",fertid)*Fert_C*0.5;
//        Cinput[hum]= Cinput[hum]+par_rothc_fert_pool_distri.get_val("HUM",fertid)*Fert_C;
//        return 0;
// }
//
//int RothC::Fill_Cinput_with_EWR(int cropid, double EWR)
//  {
//        Cinput[dpm]= Cinput[dpm]+ RothC::par_rothc_crop_pool_distri.get_val("DPM",cropid)*EWR;
//        Cinput[rpm1]=Cinput[rpm1]+RothC::par_rothc_crop_pool_distri.get_val("RPM",cropid)*EWR*0.5;
//        Cinput[rpm2]=Cinput[rpm2]+RothC::par_rothc_crop_pool_distri.get_val("RPM",cropid)*EWR*0.5;
//        Cinput[hum]= Cinput[hum]+ RothC::par_rothc_crop_pool_distri.get_val("HUM",cropid)*EWR;
//        return 0;
// }